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Abstract 

We examine the reliability of the merger trees generated for the Monte-Carlo modeling of galaxy 
formation. In particular we focus on the cold gas fraction predicted from merger trees with different 
assumptions on the progenitor distribution function, the timcstep, and the mass resolution. We show that 
the cold gas fraction is sensitive to the accuracy of the merger trees at small-mass scales of progenitors 
at high redshifts. One can reproduce the Press-Schechter prediction to a reasonable degree by adopting a 
fairly large number of redshift bins, N stC p ~ 1000, in generating merger trees, which is a factor of ten larger 
than the canonical value used in previous literature. 
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1. Introduction 

Understanding the formation and evolution of galax- 
ies is a fundamental step in linking the initial condition 
of the universe and the cosmological observational data. 
Recent systematic studies of high-redshift objects, such 
as quasars and Lyman-brcak galaxies, should provide im- 
portant clues to the early universe, although their proper 
interpretation is often not so straightforward, mainly be- 
cause those objects certainly do evolve in time. 

A theoretical study of galaxy evolution, especially its 
spectroscopic evolution, from a cosmological context, was 
begun by Tinsley (1980) and followed by many authors 
(e.g., Bruzual 1983; Arimoto, Yoshii 1986; Guidcrdoni, 
Rocca-Volmcrange 1987; Chariot, Bruzual 1991; Bruzual, 
Chariot 1993; Kodama, Arimoto 1997). These studies are 
based on a so-called 'one-zone' model which assumes that 
a galaxy does not interact with other galaxies. It is now 
fairly established, however, that structures in the universe 
have built up hierarchically from small to large scales as 
in a cold dark matter (CDM) model. This means that a 
galaxy interacts and sometimes merges with other galaxies 
even if it was an isolated system at birth. The predictions 
in the one-zone model therefore may be significantly dif- 
ferent from what happened to galaxies in a hierarchical 
universe. 

White and Frenk (1991) developed a detailed analytic 
formalism to describe the formation and evolution of 
galaxies while taking account of the hierarchical merg- 
ing of dark-matter halos, gas cooling, star formation, and 
supernova feedback. Subsequent numerical approaches in 
modeling hierarchical merging of dark halos employ two 
somewhat different algorithms; one is called the 'block 
model' in which a random-Gaussian density fluctuation 



field is generated by dividing a hypothetical rectangular 
box recursively (Cole, Kaiser 1988; Cole 1991; Cole et 
al. 1994). While this algorithm is simple and straight- 
forward, the resulting halo masses are necessarily binned 
in discrete steps of a factor of two. The other generates 
a realization of halo merger trees according to a prob- 
ability distribution function predicted by the extended 
Press-Schechter theory (Bower 1991; Bond ct al. 1991; 
Kauffmann, White 1993; Somerville, Kolatt 1999; Sheth, 
Lcmson 1999). The latter is widely used in studying 
the cosmological evolution of galaxies in a hierarchical 
universe (Kauffmann et al. 1993; Baugh et al. 1998; 
Somerville, Primack 1999; Cole et al. 2000; Nagashima, 
Gouda 2001). Throughout the present paper, we call the 
latter method the Monte-Carlo modeling of merger histo- 
ries (simply, the Monte-Carlo modeling), while it is usu- 
ally referred to as a semi-analytic model of galaxy forma- 
tion (SAM). 

The most important ingredient in Monte-Carlo model- 
ing is the conditional joint-probability distribution func- 
tion of a set of progenitor halos of mass M% at a redshift 
of Z2, which is a part of a parent halo of mass M\ at z\, 
conceptually written as 

Prob( M I , M\ ,---,M£ f ,z 2 \M 1 ,z 1 ) dM\ dM\ ■ ■ ■ dM? 

(i\T = l,.-,oo). (1) 

Unfortunately only an analytical expression for the 
conditional one-point probability distribution function, 
Prob(M2, z-i\M\, zi), is known based on the extended 
Press-Schechter theory (for the special case of the Poisson 
initial power spectra, see a different approach by Sheth, 
Lcmson 1999); one thus needs to employ an additional as- 
sumption in generating realizations of merger trees of ha- 
los in general (e.g., Kauffmann, White 1993; Somerville, 
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Kolatt 1999). Furthermore, any numerical procedure to 
generate them necessarily involves several ad hoc parame- 
ters due to the limitation of the available computation re- 
sources including the finite timestep of computation, the 
minimum mass of halos to be included in merger trees, 
and the maximum number of progenitors for each halo at 
each step. 

The purpose of the paper is to perform a systematic 
investigation of possible artificial effects of the above- 
mentioned problems on merger tree realizations, and to 
re-examine the validity of the Monte-Carlo modeling. In 
particular, we focus on the extent to which the resulting 
merger trees reproduce the conditional one-point prob- 
ability distribution function predicted by the extended 
Prcss-Schechter theory, which directly changes the frac- 
tion of cold gas. Exactly for this reason, we adopt a con- 
ventional ACDM model with the cosmological parameters 
O = 0.3, A = 0.7, h = 0.7, ct 8 = 1.0, and Q B = 0.015/t" 2 
(e.g., Kitayama, Suto 1997; Kitayama et al. 1998), and 
neglect star formation and a feedback effect for definite- 
ness. 

2. Merger Trees of Dark Matter Halos 

2.1. Constructing Merger Trees of Dark-Matter Halos 

Our model of merging histories of dark-matter halos is 
mainly based on that of Somerville and Kolatt (1999), 
which we adopt as our fiducial choice and slightly modify 
their original scheme as follows. We begin with a halo 
of mass of M\ = Af root at a redshift of z\ = z m im and 
consider its progenitors at a slightly earlier redshift of z 2 = 
zi + Az(zi). Since the joint conditional probability for the 
progenitors [equation (1)] is not known, we choose the i- 
th progenitor halo of mass M 2 according to the one-point 
conditional probability, Prob(M2, z 2 \Mi, zi), as long as 
> Af res and the total mass satisfies 



N 



i=l 



M 2 < Mi — AAf acc (< M Tes ), 



where 



AA/ acc (< M r , 



M re8 



dM 2 M 2 



dN 
dM~ 2 



(2) 



{M 2 ,z 2 \M 1 ,z 1 )(3) 



is the expectation value of the total mass of ha- 
los smaller than the resolution mass (Af ros ) with 
dN/ dM2(M2,Z2\Mi, zi) being the appropriate conditional 
mass function [equation (5) below]. In other words, we dis- 
tinguish the discrete merging and the continual accretion 
at mass M rcs , and do not resolve the halos below M res 
in our merger trees. Once all relevant progenitor halos 
are selected, we repeat the above procedure recursively 
for each progenitor until the maximum redshift (z m ax)- 
Unless otherwise stated, we set z m - ul = and z max = 15 
in the present paper. For convenience, we list in table 1 
variables which are extensively discussed in the present 
paper. 

In the original method by Somerville and Kolatt (1999), 
one stops selecting progenitors when M\ — ^ i=l M 2 be- 



comes less than M res , but without imposing the condition 
M 2 > Af ros . They carefully tuned the timesteps depending 
on M\ so that the resulting progenitor mass function be- 
comes close to equation (5) below. Rather, we stop choos- 
ing the progenitor when Mi — AM acc (< M res ) — Y2i=i M 2 
becomes negative, and the last selected progenitor M 2 
is not included in the tree. In this case, the remaining 
mass Mi — AM acc (< A/ rcs ) — Y^i=i M 2 is not necessarily 
smaller than M ros . We find that our method reproduces 
equation (5) even with the Afi-independent timesteps. 

2.2. Conditional Probability Distribution Function 

The most important and subtle issue is the proper 
choice of the one-point conditional probability, Prob(Af2, 
z 2 \Mi, zi). Bower (1991) and Bond et al. (1991) derived 
the conditional probability of M 2 at z 2 , which is a part of 
halo Mi at z\\ 



dP 
dM- 



(M 2 ,z 2 \Mi,zi) 



'Jc,2 



»c.l 



y/2n{S 2 - Si) 3 

(S c , 2 - <5 c ,i) 



exp 



2(S 2 -Si) 



dS 2 



dM 2 



(4) 



where S c j ~ 3(127r) 2 / 3 /20D(zi) (its useful approximate 
formula may be found in Kitayama, Suto 1996) is the 
critical over-density of the mass density field at a redshift 
of Zi, D{zi) is the linear growth rate, and Si = a 2 (Mi) is a 
mass variance of the density field top-hat smoothed over 
the mass scale Mi. Since equation (4) is the mass-weighted 
probability for M 2 , it is easily translated to the number- 
weighted probability that we need in the halo number 
counting: 



dN 
dM- 



(M 2 ,z 2 \Mi,zi) 



Mi dP 



(M 2 ,z 2 \Mi,zi) 



(5) 



M 2 dM 2 

Figures 1 and 2 present how the progenitor number dis- 
tribution of the merger tree realizations reproduces the 
theoretical prediction: Ajf root = 1.3 x 10 11 A/ Q (Left) and 
1.3 x 10 14 Mq (Right). In these plots, we adopt the loga- 
rithmically equal timestep in redshift : 



z (i) — (1 + ^min) X 



1 



1 



■,Wstep), 



(6) 

where iV s tep is the total number of the redshift bins. We 
defer the discussion concerning the choice of N step to the 
next subsection (2.3), and fix iV st0 p = 100 throughout this 
subsection. 

The top panels in figure 1 show the result for the 1st 
timestep (i = 1) corresponding to z {l) = 0.028 accord- 
ing to equation (6). We call this model SK-n indicat- 
ing the Somerville and Kolatt (1999) method with the 
number-weighted probability. The symbols indicate the 
average (M 2 /N ens )(AN 2 / 'AM 2 ) with the quoted error- 
bars being the corresponding onc-sigma dispersion, where 
AiV 2 is the number of progenitors in the range of mass 
M 2 - M 2 + AM 2 , and we adopt Alog 10 M 2 = 0.1. In the 
SK-n model we generate the random numbers according 
to the number-weighted probability distribution function 
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Fig. 1. Progenitor number distribution at the 1st timcstcp 
[z(i) = 0.028 with TVstep = 100 in equation (6)] for 
A/ r oot = 1-3 X 10 11 Mq (Left) and 1.3 X 10 14 M s (Right). 
Top (SK-n) and middle (SK-m) panels adopt the num- 
ber-weighted and mass-weighted conditional progenitor prob- 
ability functions [equations (5) and (4)] according to 
Somcrville and Kolatt (1999), while bottom (SL-m) panels use 
the mass- weighted one according to Sheth and Lemson (1999). 
The theoretical prediction in the extended Press-Schcchtcr 
theory [equation (5)] is plotted in solid curves. The solid tri- 
angles and open circles indicate the averages over N enB = 10 4 
merging tree realizations with the different mass conser- 
vation prescriptions, equations (2) and (7), respectively. 
The quoted error bars represent the Poisson error in each 
mass bin, and the arrows indicate the value of M TeB . 

[equation (5)]. Also we have to set the lower limit on 
the progenitor mass in adopting the SK-n model so as to 
avoid a divergent total probability. We adopted the lower 
limit of 10 _3 M rcs , and made sure that the mass function 
of progenitors of the mass range of our interest M > M ros 
is properly reproduced. 

The solid triangles show our result based on the al- 
gorithm outlined in the previous subsection. Somewhat 
surprisingly, they are completely different from the the- 
oretical distribution that we use in generating the trees 
(solid curve). Note that figure 1 plots the number distri- 
bution multiplied by M 2 , M 2 dN/dM 2 = M x dP/dM 2 [sec 
equation (5)]. To understand the origin of the discrep- 
ancy, we generate the progenitors at the 1st timestep for 
-/V cns realizations simultaneously as long as they satisfy 

TV' 

M2 < N C ns [Mi - AM acc (< M rcs )] , (7) 

i=l 

instead of the mass conservation [equation (2)] for each 
individual parent halo. In the above, N' is not the number 
of progenitors for a single halo at z = z m - mi but for an 
ensemble of N cns halos with the same mass M root . The 
resulting distribution is plotted in open circles, and in fact 
shows good agreement with the theoretical curve. 



This is simply because we attempt to generate a joint 
distribution of progenitors with a repeated use of the con- 
ditional probability [equation (5)] incorrectly; except for 
the first progenitor, the mass conservation for each halo 
[equation (2)] introduces an additional cutoff at higher 
mass in the selection probability of progenitors. In fact the 
conditional probability [equation (5)] for (z 2 — Zi)/z\ <C 1 
is sharply peaked at a mass scale M 2 just below the parent 
mass Mi , and thus even a small value of the first progen- 
itor mass may effectively bias not to choose remaining 
progenitors in the peak. Thus, the resulting distribution 
is significantly biased toward low-mass objects, i.e., the 
number density of the low-mass objects exceeds the the- 
oretical predictions by an order of magnitude (top panels 
in figure 1). 

One way out of this problem is to generate many 
(> 100) realizations simultaneously, as Kauffmann and 
White (1993) adopted. Even in this case, one needs to 
specify an additional assumption on how to plant a set 
of progenitors in a single merger tree by hand. Moreover, 
a practical implementation of this method requires one 
to discretize the halo mass, and thus becomes computa- 
tionally demanding as both the mass and time resolutions 
increase. 

Another possibility is to artificially distort the input 
conditional probability so that the selected progenitors 
obey the distribution [equation (5)]. While the required 
correction may be a fairly definite mathematical prob- 
lem, we do not know the exact answer, and thus have 
to proceed in a phenomenological fashion. Basically this 
is the approach taken by Somerville and Kolatt (1999) 
and Sheth and Lemson (1999), who adopted the mass- 
weighted probability [equation (4)] as the theoretical in- 
put. The middle and bottom panels in figure 1 show the 
resulting distribution for SK-m (Somerville, Kolatt 1999- 
mass weighted) and SL-m (Sheth, Lemson 1999-mass 
weighted) , respectively. Clearly the resulting distributions 
(filled triangles) become much closer to equation (5) un- 
der the constraint [equation (2)] although their original 
distributions [i.e., without the constraint (2)] plotted in 
open circles are completely different. 

As this indicates, the input conditional probabil- 
ity for the current purpose should be small at lower 
mass scales of M 2 relative to equation (5). Thus, we 
also attempted to make the probability proportional to 
(Mi/M 2 ) a (dP/dM 2 ). Note that the mass- and number- 
weighted probabilities [equations (4) and (5)] correspond 
to a = and 1. We were not able to obtain a similar de- 
gree of agreement for a value of a very different from 0, 
but did not find a significant change for — 0.2<a;$+0.2. 
Thus, we decided to adopt a = (the mass- weighted prob- 
ability) as Somerville and Kolatt (1999). This choice has 
an advantage that a numerical routine to generate ran- 
dom numbers becomes easier than the cases of a 0. 
Define x 2 = {5 C ,2 — S c ,i)/VS2 — Si that obeys Gaussian 
distribution of a unit variance. Equation (4) implies 
that the desired distribution of M 2 can be simply given 
via S 2 = S(M 2 ) = Si + (5 C , 2 - 5 c> i) 2 /a;|. Since Gaussian- 
distributed random numbers can be implemented easily, 
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all we have to do is to supply the inverse function of the 
mass variance, M% = S~ 1 (S2). Our SK-m implementation 
seems to yield a larger discrepancy between theory and the 
merger tree realizations at small M2 regimes than their 
original results. This may be due to the different choice 
of the timestep and the condition how to stop selecting 
progenitor halos. In any case, however, this discrepancy 
rapidly fades away in constructing the merger tree using 
many timesteps, as we show in figures 2 and 4. 




10 s 10'° 10" 10" 10>° 10" 10 lz 10 13 10'" 
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Fig. 2. Same as figure 1 except at the 25th timestep 
(z(25) = 0.028) for the mass conservation pre- 
scription [equation (2)]; Af root = 1.3 X 10 11 Mq 
(Left) and 1.3 X 10 14 M (Right). 

Figure 2 plots a snapshot of the progenitor distribu- 
tion at the 25th timestep (z = z (25) = 1.0), which makes 
sure that the mass- weighted probability reasonably works, 
even when we trace the merger tree by many steps. Also 
SK-m works a bit better than SL-m especially at small 
mass scales. Originally SL-m was proposed to correct for 
the halo exclusion effect, but does not work so efficiently 
at least in the range of parameters we surveyed. Figure 2 
also indicates that SK-n is substantially different from the 
analytical solution. This is due to the fact that SK-n tends 
to select relatively less massive progenitors preferentially 
(see figure 1), and this tendency simply accumulates in 
many steps. On the contrary, the behavior of SK-m and 
SL-m becomes closer to the analytical solution than in 
the case of figure 1. This is because the latter two models 
well approximate the probability distribution around Mi, 
which is the most important range when constructing real 
merger trees with many timesteps. 

2.3. Timestep 

The next question that we address is the appropriate 
choice of the timestep. While this is an equally impor- 
tant problem in the Monte-Carlo modeling, the previous 
authors did not discuss it in an explicit manner. 
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Fig. 3. Comparison among various timescales as a function 
of 2 in ACDM. The short-dashed and long-dashed curves in- 
dicate the cosmic time at 2, and the dynamical time of dark 
halos just virialized at 2. The dotted curve corresponds to 
the timesteps in linearly equal intervals in 2, i.e., 100 bins 
between z m i n = and 2 max = 15. The three solid curves 
correspond to the timesteps in logarithmically equal inter- 
vals in 2; N B t ep = 10, 100, and 1000 from top to bottom. 

Obviously, the timestep needs to be smaller than the dy- 
namical timescale of halos just virialized at the redshift, 
£dyn,vir(-z), because they are the objects that serve as the 
initial condition for the Monte-Carlo modeling. Figure 3 
compares this timescale and our choice [equation (6)] for 
Agtcp = 10, 100, 1000. Also, we plot the cosmic time 
icosm(z) and the timestep corresponding to the linearly 
equal bin (Az = 0.15). Even this simple comparison in- 
dicates that the logarithmic time bin with iV s tep ^ 100 is 
required. 

While the important question is how small timescales 
one should resolve, it critically depends on the problem 
that one would like to address. Therefore, we rather ask 
how many timesteps we need to reproduce the progenitor 
distribution. 

To see explicitly how different values of N stcp affect 
the realizations of merger trees, we plot in figure 4 the 
progenitor distribution with iV ste p = 1000 at redshifts of 
= 0.003, Z(i 46) = 0.499, and z (250) = 1.0. A compari- 
son of figures 2 and 4 indicates that the average progen- 
itor distribution is indeed slightly better reproduced by 
Agtcp = 1000 than N stcp — 100, particularly at small mass 
scales. 

The difference between iV s t e p = 100 and 1000 is more 
clearly illustrated when we plot the cumulative mass frac- 
tion of progenitors of mass exceeding a threshold value of 
Mthre- More specifically, figure 5 compares the theoretical 
prediction, 
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Fig. 4. Progenitor number distribution functions at the 1st, 
146th and 250th steps from the N B t ep = 1000 merger trees us- 
ing the SK-m method; Z(i) = 0.003 (Upper), ,2(1413) = 0.499 
(Middle), and 2(250) = 1-0 (Lower). The left and right 
panels are for M root = 1.3 X 10 11 M Q and 1.3 X 10 14 Mq. 
The theoretical prediction in the extended Press-Schcchtcr 
theory [equation (5)] is plotted as solid curves. 




Fig. 5. Cumulative mass fraction of progenitors of mass ex- 
ceeding a given threshold, M thre , for M roo t = 1.3 X 10 11 Mq 
(Left) and 1.3 X 10 14 Mq (Bight). The dotted and dashed 
curves indicate the averages from N stcp = 100 and 1000 merger 
tree realizations with the corresponding Poisson error-bars. 

/•M root 

F th (>M thre ;z) = / dM — (M,z\M root ,z min ), (8) 
with the average from the tree realizations. 

F mo dc\(> M thlc ;z) = — — — ^ Mhaio(z), (9) 



forM root = 1.3xlO n Af Q (Left) and 1.3x 1O 14 M (Right). 
In all panels shown here, the results with A^ tcp = 1000 bet- 
ter reproduce the theoretical prediction, mainly because 
of the small-scale behavior (see figures 2 and 4). Since 
these small mass progenitors at earlier redshifts signifi- 
cantly contribute to radiative cooling and thereby subse- 
quent star formation in the entire halo, this difference is 
indeed critical in the Monte-Carlo modeling of galaxy for- 
mation. We discuss the effect on gas cooling explicitly in 
the next section. 




z 



Fig. 6. Cumulative mass fraction of progenitors with 
different values of N stcp for Afroot = 1.3 X 10 14 Mq. 

We interpret the above as an empirical result due to 
the balance between the mass-weighted probability and 
the timesteps. If the use of the mass-weighted probabil- 
ity were strictly justified, the proper realizations with a 
smaller timestep would become more difficult numerically, 
and there is no reason why we could obtain better agree- 
ment with a larger N stcp . On the other hand, we under- 
stand that the mass-weighted probability is nothing but 
a phenomenological remedy of the problem, and with this 
choice N stcp = 1000 seems to work better than N step = 100 
empirically. 

In fact, still larger values of iV st0 p do not necessarily 
improve the result. Figure 6 is a similar plot as figure 5 
for M ro ot = 1.3 x 10 14 Mq, but with increasing iV ste p. The 
cumulative mass fraction for M > 10 11 Mq is almost un- 
changed, but the contribution from smaller mass progen- 
itors steadily increases as iV s tep becomes larger. This re- 
flects the fact that the empirical use of the mass- weighted 
conditional probability does not guarantee convergence of 
the result with respect to N step . 

We thus conclude that N stcp ~ 1000 is the optimal value 
to reproduce the Press-Schechter mass function in our 
method. 
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2.4- Number of Progenitors 

Finally, we briefly discuss how many progenitors, -/V prog , 
one should keep in order to properly reproduce the merger 
trees. Figure 7 displays the distribution functions at z = 1 
for the merger tree of M root = 1.3 x 10 14 M Q at z = 0. 
In this particular example, we use iV stC p = 1000 and the 
results are averaged over N ens = 100 realizations for each 
-Wprog- Obviously a smaller value of N prog does not prop- 
erly link the merger tree back to higher rcdshifts, and 
the number of small-mass halos is systematically under- 
predicted compared with the extended Press-Schechter 
model (solid curve). While we do not set any upper limit 
on -/Vp r0 g, figure 7 indicates that N prog ^ 5 is acceptable 
given the accuracy of the present scheme. Although some 
authors employ a binary merger tree in Monte-Carlo mod- 
cling, that scheme needs to be adjusted with a careful 
choice of the timestep and other parameters. 
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Fig. 7. Effect of the upper limit on the number of pro- 
genitors, N prog , on the progenitor distribution function 
for M root = 1-3 X 10 14 A/ at 2= 1.0 from N cnB = 100 
merger tree realizations with N B t C p = 1000. The solid 
line is the number-weighted probability function, equa- 
tion (5). The filled circles indicate the results with 
all progenitors satisfying the mass conservation [equa- 
tion (2)], while crosses, filled triangles, and open cir- 



cles are for jV D 



10, 



and 2, respectively. 



In conclusion, we have found that a reasonable agree- 
ment between the theory and the merger tree realizations 
can be obtained by employing the mass-weighted condi- 
tional probability into the Somerville and Kolatt (1999) 
scheme with A^ step ~ 1000. 

3. Gas Cooling 

So far, we have restricted our discussion to the gravita- 
tional aspect of the halo evolution. We next consider how 



the resulting merger trees affect the cold gas fraction in 
an individual halo. Of course, we must eventually discuss 
the effects on the efficiency of star formation, but we fo- 
cus on gas cooling alone, since modeling star formation, 
the feedback from supcrnovae, the chemical evolution and 
so on necessarily introduce additional (numerical) param- 
eters and the interpretation becomes more complicated. 
Thus, the principal aim of this section is to see if we can 
achieve a convergence of the cold gas fraction from differ- 
ent realizations of the merger trees. 

3.1. Description of Gas Cooling 

Our prescription of gas cooling in the merger trees goes 
as follows. First, we assume that the density profile of 
dark halos obeys the universal shape (Navarro et al. 1996): 



Phaio(r;M) 



p(z)8 c 



(r < r vir ) 



(r/r s )(l + r/r s ) 2 v ^' vir/ (10) 
(r>r vir ), 

where p(z) = f2oPco(l + z ) 3 is the mean density of the uni- 
verse at z, p c o is the present critical density, S C (M) is the 
characteristic density excess, and r vrr (M) and r s (M) in- 
dicate the virial radius and the scale radius of the halo, 
respectively. 

The virial radius is defined according to the spherical 
collapse model as 

1/3 



r vir (M) = 



/ 3M 



V47rpA nl 



(11) 



and useful approximation for the critical over-density, 
A n i = A n i(f2o,Ao), may be found in Kitayama and Suto 
(1996). The two parameters, r s and r vir , are related in 
terms of the concentration parameter, 



c{M,z\ 



_ r viv (M,z) 



[12) 



r s {M,z) ' 

We use an approximate fitting function from the simula- 
tion data of Bullock et al. (2001), 

c{M , l} = J^(J^)-°' 3 . (13 , 

be equal to 



1 + z \10 14 MqJ 

The condition that the total mass inside r vlr 
M relates S c to c as 

Sc. 



— • (14) 

3 ln(l + c)-c/(l + c) v ; 

Makino et al. (1997) showed that if the hot gas is 
isothermal and in hydrostatic equilibrium, the gas density 
profile is well approximated by the isothermal /3-model, 

Phot.O 



Phot(r) 



[1 + (r/r c ) 2 ]3/3/2 ' 



(15) 



where r c ~ 0.22r s . Wc fix (3 — 2/3 for simplicity. The 
amplitude phot,o is computed so as to reproduce the total 
hot gas in the halo when integrated up to r = r Ylr . The 
hot gas is gradually converted to cold gas according to the 
prescription below, but still the total baryon (hot + cold) 
fraction within the virial radius of each halo is set to the 
cosmic average, Hb/^o- 
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Once the gas profile is specified, one can compute the 
cooling timescale at radius r from the center of the halo, 



cool — 



3 Phot (r) k B T t 
2 /im p A(T § 



(16) 



where fi is the mean molecular weight, m p the proton 
mass, fee the Boltzmann constant, A(T gas ) the radiative 
cooling function for gas of temperature Xg as , and n^{r) the 
number density of hydrogen (including both neutral and 
ionized). We assume that the gas has the primordial abun- 
dance of hydrogen and helium (X = 0.76 and Y = 0.24, and 
thus 7Jh = Phot/ Xrrip), and compute the corresponding 
cooling function A(T) (e.g., Sasaki and Takahara 1994). 
Since we neglect the molecular and metal cooling, A(T) = 
at T < 10 4 K. 




Fig. 8. Mass of halos as a function of the virial tem- 
perature at different redshifts in the ACDM model. 

We further assume that the temperature of the hot gas, 
Tgas, is equal to the virial temperature, T v j r , of the halo. 
The relation between the virial temperature and the mass 
of a halo is plotted in figure 8. We can then solve equation 
(16) for the cooling radius, r coo i, within which the gas can 
cool within a given cooling timescale (r coo i): 



^"cool (Tcool) 



22 Vyil ( M ' Z \ pMA^viiQPhot.QTcool _ 1 ^ 



c{M,z) 



3m p /c B T v 



It now remains to define the origin of r coo i. Actually, this 
is fairly arbitrary in a sense, and we adopt the following 
simple picture. When a halo of mass Mf forms at the 
formation rcdshift, zt, its hot gas is supposed to reach 
the profile [equation (15)] instantaneously. This is defined 
to be the origin of t coo i for the halo. In the subsequent 
timesteps, we neglect the change in the hot gas profile even 
if the halo mass (M) grows due to mergers and the cooling 
radius is computed with r coo i set to the elapsed cosmic 
time since Zf. When M exceeds 2Mf, the halo is replaced 
by a newly formed massive halo, and the hot gas profile 
is reset to the profile [equation (15)] corresponding to the 



new mass and the virial temperature, and we reset the 
origin of r coo i as the new formation epoch. Incidentally 
we made sure that the value of 1.5Mf instead of 2Mf does 
not change the result, which is consistent with the finding 
of Cole et al. (2000). We apply this procedure for all 
halos, and the cold gas in each progenitor halo is simply 
accumulated (without reheated) according to the merger 
trees. 

3.2. Cold Gas Fraction in the Monte-Carlo Realization 
of Merging Histories 

In a practical implementation of the merger tree algo- 
rithm, one has to stop tracing the progenitors of halos 
of mass below the resolution mass, M res . We discuss the 
relevance of our choice of M res by looking at the cold gas 
fraction. Previous authors often apply the cutoff at a 
fixed mass or a circular velocity of halos; for instance, 
Cole ct al. (2000) consider halos with M > 5 x 10 9 /i _1 M 
in their merger trees, while Somcrvillc and Primack (1999) 
take account of halos with the circular velocity exceeding 
40 km s _1 (corresponding to the virial temperature T vir ~ 
6 x 10 4 K). The latter condition comes from an estimate 
of the smallest scale of halos which can cool in the pres- 
ence of the UV background (e.g., Thoul, Weinberg 1996; 
Kitayama, Ikeuchi 2000). 

In our present analysis, the cooling function, A(T), van- 
ishes below T = 10 4 K, since we neglect both the UV heat- 
ing and the metal/molecular cooling. Thus, we set the 
resolution mass, M res (z), as M(r v ; r = 10 4 K). As figure 8 
shows, this scale increases rapidly with time, resulting in 
a significant improvement of the computing time. On the 
other hand, this might systematically underestimate the 
cold gas fraction, since halos of mass below M res (zi) may 
have progenitors of mass larger than M res (z2) at the ear- 
lier epoch (z2 > z\). Fortunately, this is not an important 
effect, as we show below. 

To see this in detail, we plot in figure 9 the cold gas 
fraction averaged over all progenitors at z of a root halo 
of mass M root at z = 0, 

/ co id(z;Af root ) = 7 ^ — J2 M coid(M prog ),(18) 



^B-Mroot 



M prog >M reB 



for a merger tree with A^ top = 1000. If we adopt a con- 
stant value for the resolution mass, M res < 10 8 M© yields 
the convergent result for the cold gas fraction. Exactly 
the same convergence is obtained for the time-dependent 
M Tes (z) when the value is set to M(T vir = 10 4 K), but not 
if we use M(T v i r = 5x 10 4 K), for instance. This critical 
value is expected to vary depending on the thermal his- 
tory of the universe, but the appropriate value for T v j r is 
straightforwardly read off from the relevant cooling func- 
tion. Actually, M res (z) increases in this case and exceeds 
10 9 Mq (see figure 8). Thus, the required merging tree is 
less demanding from a computational point of view than 
that for M res = 10 8 Mq, for instance, as illustrated in ta- 
ble 2. Thus we decide to choose M rcs (z) = M(T vil = 10 4 K) 
for gas with the primordial abundance. 

Finally, we show the convergence with respect to N s t cp . 
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Fig. 9. Cold gas fraction [equation (18)] with different 
resolution mass, M res , averaged from N cnB = 10 merg- 
ing tree realizations at N B t cp = 1000. Upper and lower 
panels show the results for M roo t = 1.3 X 10 11 Mq 
and 1.3 X 10 14 M©; M reB (z) = M(T vir = 10 4 K), 
M(T vir = 5 x 10 4 K), 10 s M Q , and 10 9 M Q . 

Figure 10 plots the cold gas fraction for M root = 1.3 x 
10 11 M Q (Upper) and 1.3 x 10 14 M Q (Lower) for merg- 
ing trees with different A^t op . The results are fairly in 
agreement for small M root , but are very different for large 
M TOOt . This is because the merger tree at small mass 
scales, especially at M(T vlr = 10 4 K) < M < 10 10 Af Q , is 
well reproduced only when we use N stcp = 1000 (see fig- 
ure 5). When we repeat the same calculation sampling 
every 10 steps from the iV s tep = 1000 tree, the result is 
almost indistinguishable. Thus, we conclude that the pro- 
genitor distribution at small scales is quite essential in the 
estimate of the cold mass fraction of large halos. 

Incidentally the use of the timestep much smaller than 
idyn,vir(z) enables one to describe the collapse and gas 
cooling more realistically than the instantaneous approx- 
imation. While we do not attempt this in the present 
paper, this would improve the estimate of the cold gas 
fraction quantitatively. 

4. Conclusions and Discussion 

We attempted several convergence tests of the merger 
trees generated with the Monte-Carlo method. While this 
method provides a useful tool for modeling galaxy forma- 
tion in a complementary maimer to more intensive cos- 
mological simulations with ad hoc recipes of galaxy for- 
mation (e.g., Cen, Ostriker 1992; Weinberg et al. 1997; 
Yoshikawa et al. 2001), the lack of an explicit expression 
for the joint distribution function of progenitors [equa- 
tion (1)] requires one to put an additional assumption in 
practice. We confirmed that a repeated use of the mass- 
weighted conditional probability [equation (5)] reasonably 



M too[ =1.3xl0» M s 
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z 

Fig. 10. Cold gas fraction averaged from N enB = 10 merger 
tree realizations with different N B t cp . The upper and 
lower panels show the results for M root = 1.3 x 10 11 Mq 
and 1.3 x 10 14 Mq. The resolution mass of the 
tree is fixed as M rcs (z) = M(T vir = 10 4 K). 

reproduces the progenitor distribution predicted in the ex- 
tended Press-Schechter theory if one adopts fairly small 
timesteps in redshift, N step ~ 1000, a factor of ten larger 
than a typical value used in previous work. We note, 
however, that one can alternatively achieve a similar re- 
sult by fine-tuning the timestep as a function of M\ (e.g., 
Somerville, Kolatt 1999) instead of equation (6), as we 
adopted here. 

One may avoid the above problem also by using merger 
trees generated via iV-body simulations (Kauffmann et 
al. 1999a; Somerville et al. 2001). In fact, they claim that 
the agreement between the A^-body simulations and the 
Monte-Carlo method is good. Benson et al. (2001) com- 
pared the SPH simulations and the Monte-Carlo model- 
ing, and concluded that both agree with each other on 
the cold gas mass fraction and mass function of the ha- 
los. While this comparison is encouraging, it is not yet 
clear if the lack of the joint distribution function of pro- 
genitors [equation (1)] in the Monte-Carlo modeling may 
not be essential. Thus, further detailed studies are def- 
initely important to test the reliability of both TV-body 
and the Monte-Carlo modeling in generating merger tree 
realizations. 
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Table 1. Summary of the variables used in building merger trees and in gas cooling. 



Symbol 


Adopted value 


Physical meaning 


M root 




mass of halo at z = z m in 


T ■ 

1 vir 




virial temperature of halo 


M rcs 


M(T vir = 10 4 K) 


minimal mass of progenitors resolved in each merger tree 


7"cool 


halo mass doubling time 


cooling time scale for gas in hosting halos 




1000 


number of redshift bins (logarithmically equal interval) 


N ons 




number of realizations of merger trees 


-2-min 





minimum redshift of merger trees 


^max 


15 


maximum redshift of merger trees 



Table 2. CPU timing of the Monte-Carlo modeling for one merger tree on a 21264 alpha 600 MHz machine. 



-™root 


N st cp 




Number of progenitors 


CPU-time ( 




(M e ) 








merger tree gas 


cooling 


1.3 x 10 11 


100 


10 8 Mq 


3876 


3.4 


0.55 


1.3 x 10 11 


100 


M(T vir = 10 4 K) 


1687 


2.9 


0.43 


1.3 x 10 11 


1000 


1O 8 M 


56057 


15.1 


4.4 


1.3 x 10 11 


1000 


Af(T vir = 10 4 K) 


20556 


6.8 


2.3 


1.3 x 10 i4 


100 


1O 8 M 


1122420 


244 


75.6 


1.3 x 10 14 


100 


M(T vir = 10 4 K) 


640409 


137 


46.5 


1.3 x 10 14 


1000 


10 8 Mq 


29725672 


6514 


1993 


1.3 x 10 14 


1000 


M(T vir = 10 4 K) 


18875982 


4140 


1400 



